We present an original algorithm to simulate seismic wave propagation in anisotropic elastic media. The algorithm is essentially oriented towards solving real-life applications. It is designed to simulate multi-shot multi-receiver seismic data for anisotropic seismic models. Implementation of the algorithm is based on three parallelization levels: direct parallelism to the right-hand sides, domain decomposition-based parallelism governed under MPI functionality, and stencil computations which are easy to perform with multi-core GPU architecture. We illustrate its applicability and high performance in solving real-life problems by simulating two full-coverage datasets consisting of 6600 sources each.